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Abstract. We simulate the neutrino and 7-ray emissions of the Galaxy which are 
originated from the hadronic scattering of cosmic rays (CR) with the interstellar 
medium (ISM). Rather than assuming a uniform CR density, we estimate the spatial 
distribution of CR nuclei by means of numerical simulations. We consider several 
models of the galactic magnetic field and of the ISM distribution finding only a weak 
dependence of our results on their choice. We found that by extrapolating the predicted 
7-ray spectra down to few GeV we get a good agreement with EGRET measurements. 
Then, we can reliably compare our predictions with available observations above the 
TeV both for the 7-rays and the neutrinos. We confirm that the excesses observed by 
MILAGRO in the Cygnus region and by HESS in the Galactic Centre Ridge cannot be 
explained without invoking significant CR over-densities in those regions. Finally, we 
discuss the perspectives that a km 3 neutrino telescope based in the North hemisphere 
has to measure the diffuse emission from the inner Galaxy. 



1. Introduction 

It is well known that the interaction of cosmic ray (CR) nuclei with the interstellar 
medium (ISM) should give rise to diffuse 7-ray and neutrino emissions [H [2] . Several 
satellites (see [3 J for a review), especially EGRET [3J [5], already observed the 7-ray 
diffuse emission of the Galaxy up to ~ 10 GeV. Around the GeV, the hadronic origin of 
a large fraction of that emission is testified by its good correlation with the distribution 
of the interstellar hydrogen. GLAST observatory [6] should soon provide a deeper 
insight into the nature of that emission as it will probe the 7-ray sky with much better 
sensitivity and angular resolution than EGRET and extend the explored energy window 
up to 300 GeV. Above that energy, however, satellite observatories can hardly provide 
significant data. 

Ground based experiments are divided in two main groups: extensive air shower 
arrays and atmospheric Cherenkov telescopes. In the former class, experiments like 
MILAGRO [7\ and TIBET [8] can probe large regions of the sky at energies larger than 
few TeV's with an angular resolution of several degrees. Those experiments already 
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provided interesting constraints on the diffuse emission from the galactic disk [H] HO] 
and, in the case of MILAGRO, a diffuse emission from the Cygnus region has been 
recently observed [TT]. Unfortunately, none of those experiments covers a sky window 
encompassing the Galactic Centre (GC). Atmospheric Cherenkov telescopes, like HESS 
[T2] and MAGIC [13], are sensitive to 7-rays with energy in the range 0.1 < E < 100 TeV 
and have a very good angular resolution (better than 0.1°). However, although these 
instruments can search for a diffuse emission in some limited regions of the sky, they 
are mainly dedicated to observe quite localised sources. 

Neutrino telescopes (NTs) may soon provide a new valuable probe of high energy 
CR physics which can be complementary to 7-ray observatories. The large detection 
volume, large field of view, available observation time, and relatively good angular 
resolution (which should be better than 1° for water Cherenkov NTs looking for up- 
going muon neutrinos) make, at least in principle, these instruments well suited to 
search for the diffuse neutrino emission of the Galaxy. 

The AMANDA telescope, which is operating at the South Pole, put already an upper 
limit on the u^s flux from the GP [14J (see section HD and ICECUBE [15] should 
significantly improve that limit. 

However, in this context the most promising NTs are those based in the North 
hemisphere (ANTARES [IB], NEMO p2] and NESTOR [H] and especially their planned 
km 3 upgrade, the Km3NeT project [19]) as they will be sensitive to up-going z/^'s coming 
from the GC region. Furthermore, since these instruments will be deployed in water, 
their angular resolution should be better than NT at the South pole improving their 
chances to disentangle the signal from the background. 

From the theoretical point of view, several calculations of the 7-ray and neutrino 
diffuse fluxes above the TeV have been already performed. In [201 (2TJ [22] the CR 
distribution in the Galactic disk was assumed to be homogeneous so that the secondary 
v and 7-ray fluxes come out to be proportional to the gas column density (in [20], [22] a 
uniform gas distribution was assumed). While that was a very reasonable approximation 
to start with, a more accurate analysis should account for the inhomogeneous CR 
distribution which arises due to the not uniform scatter of sources. It should be noted 
that, since CR sources are likely to be related to supernova remnants (SNR) and those 
objects are most abundant in the gas (target) rich regions, the secondary 7-ray and 
neutrino emissions may be significantly enhanced with respect to the case in which a 
uniform CR distribution is assumed. The CR homogeneity assumption was released in 
more recent works where the high energy nuclei distribution was modelled by solving 
the spatial diffusion equation [23], [24]. Those papers, however, were mainly addressed 
to model the diffuse 7-ray emission below the TeV and did not consider neutrinos. In 
[25] a different, and promising, approach was considered to solve the diffusion equation 
and to model the neutrino emission up to the PeV. In that paper, however, several 
approximations were done which did not allow to reach an angular accuracy better than 
10°. 

The aim of this work is to model both the diffuse 7-ray and neutrino emissions of 
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the Galaxy above the TeV with a better accuracy. That is necessary in order to be able 
to establish the detection perspectives of air shower arrays and neutrino telescopes and 
to interpret correctly forthcoming observations. Our approach is comprehensive as we 
carefully account both for the spatial distribution of CR nuclei and for that of the ISM. 

In section [2] we start by discussing the main properties of the ISM. We pay special 
attention to the distribution of SNR (see section 12.11) and to that of the atomic and 
molecular hydrogen (section [2.31) . In both cases we adopt models which we apply for the 
first time in this context. Details of gas models and their comparison with those adopted 
in previous works are given in the Appendix A In section[3]we model the distribution of 
the CR nuclei in the Galaxy by solving numerically the diffusion equation. Our approach 
is similar to that sketched in [25] that we fully exploit here. Differently from [231 121] we 
use expressions for the diffusion coefficients as determined from MonteCarlo simulations 
of charged particle propagation in turbulent magnetic fields [26l [27]. This allows us 
to get more detailed CR distributions above the TeV and to test how much the large 
uncertainties in the knowledge of the turbulent component of the GMF (see section 12.21) 
affect our predictions. In section H] we combine the simulated CR distributions and the 
gas models to map the expected and 7-ray emissions. Then, in section [5] we compare 
our predictions with available experimental results. In section [6] we briefly discuss the 
perspectives that a km 3 NT to be built in the Mediterranean sea has to detect the 
muon neutrino emission from the inner Galaxy. Finally, in section [7J we summarise our 
conclusions. 



2. The spatial structure of the ISM 

In order to assess the problem of the propagation of CRs and their interaction with the 
ISM we need the knowledge of three basic physical inputs, namely: 

(i) the distribution of SuperNova Remnants (SNR) which we assume to trace that of 
CR sources; 

(ii) the properties of the Galactic Magnetic Field (GMF) in which the propagation 
occurs; 

(iii) the distribution of the diffuse gas providing the target for the production of 7-rays 
and neutrinos through hadronic interactions. 

Since all these inputs are affected by large uncertainties and systematics, we need to 
estimate how their poor knowledge is transferred in the resulting 7-ray and v emission. 
Therefore, we compare the results obtained using several models. In the following we 
assume cylindrical symmetry and adopt the Sun galactocentric distance r Q = 8.5 kpc. 
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2.1. The SNR distribution in Galaxy 



The rate of galactic SN explosions, as inferred from observations of SN in external 
galaxies similar to the Milky Way |28j, is roughly 

ft ~ ^ yr- «„ ~ 1 yr- , (1) 

respectively for type I and type II SN. The total rate is ~ 1/48 yr 1 which is in 
reasonably good agreement with records of historical SN. Since we will normalise the 
CR injection rate by requiring the simulated CR flux to meet the observed value at 
Earth position, we are only interested in the R\ and Ru ratio. 

What is even more crucial to our analysis is the spatial distribution of type I and 
II SN. Since type-la and core-collapsed SN originate from different stellar populations 
their spatial distribution, as well as their rate, are different. 

Several methods to determined the SNR distribution in the Galaxy are discussed in 
the literature. One of the most commonly adopted methods is that to estimate the SNR 
distances on the basis of the surface brightness - distance (E — D) relation [29]. Such 
an analysis does not cover the GC region. Furthermore, several doubts have been risen 
on the accuracy of that method as it is plagued by a number of systematics concerning 
the completeness of the available SNR catalogue and the proper handling of selection 
effects [30l [3T] . Here we adopt a SNR distribution a distribution which is inferred from 
observations of related objects, such as pulsars or progenitor stars, as done e.g. in [32] 
which, in our opinion, is a safer approach. 

In [32] the spatial distribution of type-la SNR adopted was assumed to follow that 
of old disk stars which have an exponential scale length ~ 4.5 kpc along r and an 
exponential scale height ~ 300 pc. Therefore 

r — r & \z\ \ 



TZ\(r } z) = Kj exp 



(2) 



4.5 kpc 0.3 kpc^ 

where K\ is a normalisation factor. Although type-la SN are globally less frequent than 
core-collapsed SN, their rate is dominating in the inner few kpc's of the Galaxy. 

To trace core-collapsed originated SNR one may use either HII regions, which are 
produced by their luminous progenitors, or pulsars, which are a likely left-over of the 

EE 



collapse. The pulsar distribution at birth was estimated to be 

3.55 exp 



K n (r,z) =K n f(z) 



r — 3.7 kpc 
2.1 kpc 



exp 



r — 77 



(6.8 kpc)' 



r < 3.7 kpc 



r > 3.7 kpc 



(3) 



where 



f(z) = 0.79 exp 



0.21 exp 



0.6 kpc 



(4) 



\0.2 kpc J 

While the absolute values of K\ and K\\ are irrelevant here, their ratio is needed to 
normalise the relative weights of type-I and type-II SNR distributions. By requiring 
that Ru/Ri = 4.2, as it follows from ([I]), we find K n /Ki ~ 7.3. 
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In figure [T] we show the total SNR radial distribution as obtained in that way and 
compare it with that given in |29j. Both distributions have been normalised so to take 
the same value at the solar circle where observations are most reliable. It is evident 
that in the inner Galaxy the SNR rate distribution which we adopt in this work, is 
significantly higher than that given in [29]. Rather, we have a relative good agreement 
with the distribution adopted in [35] (see also [36]) which is based on pulsars, and with 
independent observations of the 1809 keV line of 26 Al which is thought to be a reliable 
tracer of SNs [3ZIEH]. 



4 F 




Figure 1. The SNR radial distribution f(r, z — 0) is shown in arbitrary units versus 
the distance from the GC. The distribution is such that, in cylindrical symmetry, 
2ir J rf(r, z) dr dz represents the total SN rate in the whole galaxy. The upper line 
(red, continue) is that derived in [32] which we adopt here, while the lower one (blue, 
dashed) is from [29]. They are both normalised to unit at r = r = 8.5 kpc. 

We assume the spatial distribution of the CR sources to coincide with that of SNRs 
and the CR energy spectrum at source to be everywhere a power-law with exponent (3 
(P, which refers to the injection spectrum, should not be confused with the slope of the 
observed CR spectrum). Therefore, the source term in our diffusion equation will be 

Q{E, r, z) = q{r, z) E^ = K (TZ^r, z) + TZ n (r, z)) E' 13 . (5) 

Its absolute normalisation K and spectral slope (3 will be fixed so to reproduce the 
observed spectrum at Earth below the knee. For protons [39] 

I P (E) ~ 8.7 x 1(T 6 (p^y) ( TeV cm2 s sr ) _1 • ( 6 ) 
2.2. Regular and random magnetic fields 

The Milky Way, as well as other spiral galaxies, is known to be permeated by large-scale, 
so called regular, magnetic fields. The orientation and strength of those fields is measured 
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mainly by means of Faraday Rotation Measurements (RMs) of polarised radio sources. 
From those observations it is known that the regular field in the disk of the Galaxy is 
prevalently oriented along the disk plane and it seems to follow the galactic arms as 
observed in other spiral galaxies. According to [40] , its strength at the Sun position is 

B = 5 r d «!f (r ,O) = 2.1 ± 0.3 fiG while at smaller radii B?^(r) = B exp j- 7 " ~ r ° j 

where tb = 8.5 ±4.7 kpc. A 1/r profile seems to give a worst fit of data. Unfortunately, 
observations are not significant for r < 3 kpc. Most likely [ID] the regular field in the 
disk has a bi-symmetric structure (BBS) with a counterclockwise field in the spiral arms 
and clockwise in the interarm regions. Concerning its vertical behaviour, it is generally 
assumed that B^g k decreases exponentially for increasing values of \z\ with a scale height 
of few hundred parsecs. There are increasing evidences that the field is symmetric for 
z -> -z (BBS-S) @D. 

Superimposed to the regular field a random, or turbulent, component of the GMF 
is known to be present. In the disk, this component is comparable to, or even larger 
than, the regular one. Indeed, the locally observed rms value of the total field is about 
6 ± 2 fiG. which is 2 — 4 times larger than B reg (r Q , 0). From polarimetric measurements 
of stellar light and RMs of close pulsars it has been inferred that the GMF is chaotic on 
all scales below i max ~ 100 pc. The power spectrum of the GMF fluctuations is poorly 
known. Observational data, obtained from RM of pairs of close pulsars, are compatible 
with a Kolmogorov spectrum, i.e. B 2 (k) oc k~ 5 ^ 3 , though with a very large uncertainty 
(see e.g. [12] and ref.s therein). 

What is most relevant here is the galactic magnetic halo (MH) since most of CR 
propagation takes place outside the disk. The regular component of the MH has been 
studied by means of RMs of polarised extra-galactic radio sources. Those observations 
showed that Bj^ is mainly azimutal, the same as B d g g k , and that its vertical scale 
height is z r ~ 1.5 kpc [13]. The radial behaviour of B^ is poorly known and it is 
generally assumed that it traces that of B^r ■ Here we share the same attitude. It is 
worth noticing that the vertical symmetry of the MH might be opposite with respect to 
that in the disk: this fact may have relevant consequences for the propagation of CRs 
with E/Z > 10 15 eV [27]. 

In the following we assume a symmetric structure so that the regular component 
of the MH can be considered as an extension of B d '?. k and we can combine both in a 
unique simple structure: 

B Ies (r,z) = 5 exp {-^} 2c J [z/Zr) . (?) 

with z r = 1.5 kpc and tb left as a free parameter. Similarly to [27J we use a cosh function 
to regularise the exponential vertical profile of the regular and turbulent components of 
the GMF at z = 0. 

Concerning the properties of the random MF component in the halo, very little is 
known from observations. Both the vertical extension of the radio halo and the isotopic 
relative abundances of CR suggest that the vertical scale height of the random fields 
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is significantly larger than that of the regular one, z t — 3 — 5 kpc or larger. It is 
also likely that the turbulence strength increases in the regions with the highest star- 
forming activity. We account for this possible radial dependence of -B tU rb by means of 

the parameter a(r) = - — ""^ r ™ s ( ' — 1_ Therefore, we write 

B reg (r,0) 

iWr,*) = »(r)iWr, 0)5-^. (8) 

Due to the relative small number of extra-galactic radio sources, the maximal scale 
of magnetic fluctuations L^}^ and their power spectrum in the MH are practically 
unknown. In the following; we assume that L[^° is the same as in the disk and consider 
only Kolmogorov and Kraichnan (^B 2 (k) oc k~ 3 ^ 2 ^j power spectra. 

2.3. The gas distribution 

The diffuse gas accounts for about 10-15% of the total mass of the Galactic disk (not 
including dark matter) and its chemical composition is dominated by hydrogen (about 
90.8 % by number and 70.4 % by mass) and helium (9.1 % by number and 28.1 % 
by mass). Hydrogen is shared in three main components [32]: ionised (HII, total 
mass M H ii ~ 1 x 10 9 M©), atomic (HI, M m ~ 6 x 10 9 M ) and molecular (H 2 , 
Mh 2 — 1 — 2 x 10 9 Mq). While all these components have comparable masses, their 
spatial distributions are quite different. The much hotter HII has a scale height along 
the vertical axis considerably larger than the other hydrogen components (/ihii — 1 kpc). 
Therefore its contribution to the 7-ray and neutrino emissivity from the galactic plane 
is subdominant and can be neglected in the following. 

The construction of galacto-centric radial density profiles for the HI and H 2 



intrinsically needs model dependent assumptions (see Appendix A for further details). 
Therefore large uncertainties are involved in this operation. Here we consider two 
models, which we call gas models A and B, both for the atomic and molecular hydrogen, 



that will be discussed in details in Appendix A 



Model A It has been developed by Nakanishi and Sofue (NS) in [44J for the HI and 
by the same authors in [45] for the H 2 . 

Model B We construct model B by suitably combining the results of different analyses 
which have been separately performed for the disk and the galactic bulge. For the 
H 2 and HI distributions in the bulge we use a detailed 3D model recently developed 
by Ferriere et al. [46] on the basis of several observations. For the molecular 
hydrogen in the disk we use the well known Bronfman's et al. model [47]. For 
the HI distribution in the disk, we adopt Wolfire et al. [IB] 2-dimensional model. 
Although Ferriere et al. model is not cylindrically symmetric, we verified that by 
fitting their 3D distribution with a cylindrically symmetric one and assuming a 
Gaussian vertical profile peaked on the GP, we get 7-ray and neutrino fluxes that, 
when integrated over windows larger than 1 degrees squared, differ very little from 
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those obtained using the complete 3D models. For this reason, in the following we 
work only with averaged 2D distributions. 

In figure [2] and figure [3] we show the volume density radial and vertical profiles of 
the HI and the H2 obtained with models A and B together with the continuos fits we use 
in our analysis without reporting observational errors as they are typically much smaller 
than systematic ones (the difference between model A and B will provide a glimpse of 
the amount of those uncertainties). 




Figure 2. Left panel: the HI volume density along the Galactic Plane (z = 0) is 
plotted versus the radius r. Green crosses display the binned distribution as given in 
[44] while and the continuos line (green) is our fit of those points (model A) . Blue 
crosses are derived from [46]. The dashed (red) line is our fit of those points with the 
distribution given in [48]. for r > 3 kpc (model B). Right panel: the radial profile of 
the HWHM scale height (Z1/2) for the same models. 

While the gas densities in model A and B differ relatively little close to the solar 
circle (what is most relevant here is the H2) the main discrepancies arise in the galactic 
bulge. Indeed, that is the region where the uncertainties on the gas velocity are the 
largest. That discrepancy, however, has little consequences on the gas column density 
(see figure Hj) as it should since this quantity is almost directly related to the observed 
CO emission. 

In the following we will use model B as our reference model since, in the central 
region of the GE, it provides a better fit of the 12 CO emission survey [49] . In figure [4] we 
also compare the hydrogen column density distributions obtained with models A and 
B with that used in [21] finding that the former are more narrowly peaked along the 
GP. In the following we will assume that helium is distributed in the same way as all 
hydrogen nuclei. 
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Figure 3. Left panel: the H2 volume density along the Galactic Plane as a function 
of r. Green crosses display the binned distribution as given in [45j and the continuos 
(green) line is our continuos fit of those points (model A) . Red and blue crosses are 
derived (see text) from [47] and [46] respectively. The dashed (red) line is our combined 
fit of all those points (model B). Right panel: the radial profile of the HWHM scale 
height (zi/2) for the same models. 




-50 50 -5 5 
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Figure 4. Profiles of the hydrogen nuclei (HI + 2H2) column density are shown 
along the galactic plane (left panel) and along the I = line (right panel). The three 
diagrams correspond to model A (dashed, blue line), model B (continuos, red line) 
and to the hydrogen distribution adopted in [21] . Column densities are averaged over 
1° x 1° angular bins. 



3. Numerical Simulation of CR diffusion in the Galaxy 

The ISM is a quite turbulent magneto- hydro-dynamic (MHD) environment. Since the 
Larmor radius of high energy nuclei, 

rm = - 0.1 pc , (9) 



ZeB ieg ' V10 2 TeV/ \ 1 
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is much smaller than the maximal scale of the magnetic field fluctuations L max ~ 100 pc, 
the propagation of those particles takes place in the spatial diffusion regime. The 
diffusion equation which describes such a propagation is (see e.g. [50] ) 

V^(£,r,z) = -V t (D l , i (r,z)V J N(E,r,z))=Q(E,r,z) (10) 

where N(E,r,z) is the differential CR density averaged over a scale larger than L max , 
Q(E, r, z) is the CR source term (jSJ) and Dij(E, r, z) are the spatial components of the 
diffusion tensor. We assume cylindrical symmetry and use Einstein notation for repeated 
indexes. Equation (JlO]) does not contain energy loss/gain terms. This is justified by 
the fact that in the high energy range relevant for our problem [E v ( 7 ) > 1 TeV i.e. 
(E) > 10 TeV) nuclei energy losses and re-acceleration are expected to be negligible 
(see e.g. [HUEI]). 

Under the assumption of cylindrical symmetry the diffusion equation takes the 
simpler form [50] 



{-~d r [rD ± d r )-d z [D ± d z )+u r d r + u z d z }N{E,r,z)=Q{E,r,z) , (11) 

where D± and Da are respectively the diffusion coefficient in the direction perpendicular 
to B reg and the antisymmetric (Hall) coefficient, while u r and u z are the drift velocities 



as defined in Appendix B Equation ( fiTT) will be derived in Appendix B 



Several simulations of the diffusion coefficients in high turbulence MHD media 
have been performed [521126] . We adopt here the expressions for these coefficients which 
have been derived in [53] and used in [25]. Those expressions provide D± and Da 
as functions of the regular magnetic field and of the turbulence level a in that point. 
Respect to other works, where the CR density has been simulated by using a mean value 
of the diffusion coefficients as derived from the observed secondary /primary ratio of CR 
nuclear species, our approach offers the advantage to provide the diffusion coefficients 
point-by-point. This allows to test how the expected neutrino and 7-ray emissions depend 
on the properties of the turbulent components of the GMF. We verified, by means of 
dedicated runs, that the simulated escape time of nuclei found with our numerical code is 
compatible with that estimated from observations of the B/C ratio measured at energies 
around the GeV (see e.g. [21]). 

All the advantages of getting a realistic spatial distribution of CR by solving the 
diffusion equation, however, were not fully exploited in [25]. In that work, in fact, the 
author solved equation (fill) analytically, which was possible only considering an over- 
simplified spatial dependence of the diffusion coefficients. For example, D± was assumed 
to be spatially homogeneous. 

In this work we solve numericall)||] (jTTI) under more general conditions. We impose the 
boundary conditions iV = on a the cylindrical surface (r = 30 kpc, z = z t ) where the 
turbulent halo is supposed to end and CRs escape to infinity without further diffusing. 



I The code is based on a previous Fortran version wrote by J. Candia which we translated in C++ and 
improved in several parts. 
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For any given value of the particle rigidity the code provides a bi-dimensional histogram 
i[}(r, z) mapping the particle spatial distribution in the stationary limit. 

3.1. Models 

In order to verify how much the large uncertainties in the knowledge of the GMF 
properties may affect our predictions for the CR distribution, we consider and test 
several models. 

Our first result of these tests is that the CR spatial distribution and spectrum 
are almost independent of the radial and vertical lenght-scales of the regular magnetic 
field. Therefore the large observational uncertainties in tb and z r do not affect our final 
results. In the following we will always assume re = 8.5 kpc and z r = 1.5 kpc. 

More significant is the effect of changing the turbulent halo scale height z t as that 
parameter determines the length over which CR have to diffuse before escaping to infinity 
(here we always assume that z r < z t , as suggested form observations; see section l2~2l) . 
Indeed we found that in the inner Galaxy, i.e. where sources are most abundant, the 
CR density increases with z t . We considered several values of that parameter finding 
that z t = 3 kpc provide the best matching of EGRET data (see section [5]). 

In table [T] we summarise the main features of all models that will be discussed in 
this paper. 

Table 1. The main properties of the models considered in this section. 



model # 


SNR 


o-(r) 


turbulence 


z-symmetry 





CB [29] 


1 


Kolmogorov 


S 


1 


Ferriere [32] 


1 


Kolmogorov 


S 


2 


Ferriere [32] 


like SNR 


Kolmogorov 


S 


3 


Ferriere [32] 


1 


Kraichnan 


s 


4 


Ferriere [32] 


1 


Kolmogorov 


A 



In figure [5] and [6] we show the high energy proton fluxes along two significant sections 
of the galactic halo: the galactic plane and a plane perpendicular to it at r = r . The 
effect of changing the SNR distribution is evident when comparing models and 1 
as both have the same GMF turbulent spectrum (Kolmogorov) and strength (a = 1 
everywhere on the GP). It is evident that the adoption of the SNR distribution given 
in [32] gives rise to a ~ 30 % increase in the CR density in the GC region with respect 
to that derived following [29] . 

We also investigate the effect of changing the turbulence level uniformly in the 
magnetic halo finding a marginal effect at energies below the PeV. More interesting is 
the effect of assuming a radially dependent turbulence strength cr(r). In model 2 we 
assume that quantity to follow the same radial profile of SNR as shown in the upper 
curve of figure [U That choice is justified by the well known argument according to 
which MHD turbulence in the ISM is powered by SN ejecta. By comparing proton flux 



Diffuse Neutrino and Gamma-ray Emissions of the Galaxy above the Te V 



12 



profiles in figure M for the models 1 and 2 (as both are obtained by using the same 
value of a — 1 at the solar circle) we find that a radially dependent a gives rise to a 
smoothing of the CR density distribution with respect to that of sources. That is to 
be expected as regions which are poor of sources are more easily filled by CR coming 
from more active regions if the turbulence strength is locally smaller. The effect is less 
evident at higher energies since in that case CR escape more rapidly along the z axis. 

In the above we assumed a Kolmogorov power spectrum (i.e. 7 = 5/3) for the 
turbulent component of the GMF. It is interesting to investigate how the CR spatial 
distribution is affected by adopting a different spectrum. A reasonable possibility is to 
assume a Kraichnan spectrum (model 3) which is characterised by 7 = 3/2. Noticeably, 
in this case D± oc E 1 ! 2 so that the value of the spectral index at the sources, which is 
required to fit the observed proton spectrum at the Earth, is (3 = 2.7 — 0.5 = 2.2. This is 
in good agreement with the value measured for the 7-ray spectrum of several SNRs (see 
e.g. (54] and ref.s therein). Since the Kraichnan spectrum is harder than Kolmogorov's, 
it gives rise to a more tight confinement of CR in the nearby of sources explaining the 
higher density in the inner Galaxy which we observe for that model in figure [6j 

The deformation of the vertical profiles observed for all models at high energy is 
due to the Hall diffusion which gives rise to a drift in a direction perpendicular to both 
B rog (along 0) and the radial component of V<fr p . 
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Figure 5. The proton flux profiles at E = 10 TeV are plotted along the galactic 
midplane (left panel) and along the z-axis at r = r Q (right panel). The dash-dotted 
line (orange), continuos (red), dotted (green), and dashed (blue) correspond to the 
models 0,2,1,3 respectively (see text). All fluxes have been normalised to the observed 
value at the Sun position. 



Finally, we also considered the effect of changing the symmetry of the regular 
magnetic field with respect to the galactic plane (model 4). In all previous cases a 
symmetric (S) configuration (i.e. B reg does not reverse its sign for z — > —z) was 
considered. Since observations did not settle what the actual symmetry is yet, it is 
interesting to consider the possible implications of its change. In [501 127] it was already 
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Figure 6. The same as figure [5] for E = 1 PeV. 



shown that if the magnetic halo is A-symmetric, Hall diffusion gives rise to a significant 
increase of the CR density in the GC region at energies above the PeV. We confirm the 
existence of that effect. The consequences for the neutrino spectrum where investigated 
in [25]. However, as they are negligible in the energy range considered in this work we 
will disregard the A-symmetric models in the following. 

As we mentioned, the diagrams in this sections refer to protons only. The 
corresponding flux profiles for composite nuclei can be obtained by a simple energy 
rescaling (E —> E/Z) and by using the proper normalisation given by the observed flux 
of each nuclide. Therefore, 

dn z (E; r,z) ^ , , „ , „ / E ^ "° z 



dE 



En 



^(E/Z;r } z) 



(12) 



z no- 
where fz and az are respectively the locally observed fraction at 1 TeV of the species 
with charge Z with respect to protons and their spectral slope. Recent compilations 
of measurements of both quantities can be found in [39|. We conclude this section by 
noticing that although the CR energy spectrum is not exactly the same over the whole 
galaxy, this has negligible implications for the following results. 



4. Mapping the 7-ray and neutrino emission 

We start by considering 7-ray production by the decay of neutral pions which are 
generated by the interaction of the nucleonic component of CR with the ISM (mainly 
hydrogen and helium). Here we do not consider a possible contribution to the diffuse 
7-ray flux which may be originated by Inverse-Compton (IC) scattering of relativistic 
electrons with the background radiation (see e.g. [51] [23]). Bremsstrahlung is negligible 
at the energies considered in this work. Due to the low density of the ISM practically all 
mesons decay before interacting with matter. Since at those high energies constituent 
nucleons interact independently one from the other, here we need only to consider 
elementary inelastic nucleon-nucleon scattering. Furthermore, we can safely assume 
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proton-neutron invariance. In the energy range considered in this work {E 1 < 100 TeV), 
the attenuation of the 7-ray flux due to pair-production scattering onto the CMB 
radiation can be neglected (see e.g. [21] 122]). 

High energy neutrinos (E > 1 TeV) are prevalently originated by the decay of 
charged pions and kaons which are generated by the same hadronic scattering process 
considered for the production of 7-rays. At the source only electron and muon neutrinos 
are produced. During the propagation over galactic distances neutrino oscillations 
redistribute equally the neutrino budget among all lepton families. In the following 
we will be interested only in the + flux reaching the Earth. 

Under the assumption that the primary proton spectrum is a power-law and that 
the differential cross-section follows a scaling behaviour (which is well justified at the 
energies considered in this paper), the 7-ray (muon neutrino) emissivity, can be written 
as (see e.g. [U [2]) 

Qi (^+u^){E u ; r, z) = j N dn p(E^ r ' z ) appC nH fa ^ Yl ( ^ + p M )(a) . (13) 

Here a is the primary CR spectral index and a pp ~ 3.3 x 10~ 26 cm -2 is the pp inelastic 
cross-section at 1 TeV [53]. The 7-ray (muon neutrino) yield Yy (vu+vu,)( a )i as determined 
in [56], is shown in figure ^M. Those values are in agreement with previous results (see 
e.g. [S3 ESI EH]) within 20%. Above few GeV's, yields are practically independent of 
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Figure 7. The 7-ray and + yields are represented as a function of the primary 
nuclei spectral index a. Neutrino oscillations are accounted for. 

energy. For a = 2.7 the 7/1^ + flux ratio is 3.1. 

The factor /jy in (flBl) represents the contribution from all other nuclear species 
both in the CR and the ISM. We find 




§ An almost 10% contribution to the photon emissivity coming from 77 decay, which was not considered 
in [56] . has been included here. 
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which is in good agreement with the value found in other works (see e.g. [6D| I6"T]). Here 
we used experimental values of fz and az as given in [39J and the helium/hydrogen 
ratio in the ISM n^/nn — 0.09 [32]. All other nuclear components in the ISM give a 
negligible contribution. 

The differential 7-ray (neutrino) flux reaching the Earth is given by the line integral 

J 7(^+^)(^' M) = 4~ / ^7 (vr+BriiEvl M,«) ds , (15) 

where s is the distance from the Earth and (I, b, s) are related to (r, z, <p) through 

z = ssinb r = \J (s cos b cos I — r ) 2 + (s cos b sin I) 2 , (16) 

in cylindrical symmetry. Finally, the integrated fluxes are determined by integrating 
the power law spectrum over the energy up to 1 PeV. 

For the sake of clarity, in the following we will show flux diagrams as obtained 
only with our model 3B (model 3 for the CR distribution and B for the gas) which 
is our preferred model. As we mentioned in section 12.31 the gas model B gives the 
best matching of CO surveys, while the CR model 3 has to be preferred because, by 
adopting gas model B, it best reproduces EGRET observations above few GeV's (see 
section I5TT1) . At the end of this section we will briefly discuss how our predictions would 
change adopting different models. 

In figure M we show two representative sections of the neutrino flux profile above 1 
TeV . In order to show how the expected signal may depend on the experimental angular 
resolution, in the same figure we draw the flux averaged over angular bins of different 
sizes. It is evident that due to the narrowly peaked behaviour of the gas density along 
the GP (b = 0), the averaged flux which may be measured from these regions should 
change significantly by varying the angular resolution. 




I (degrees) b (degrees) 



Figure 8. The integrated + flux is shown along the galactic equator (b = 0, left 
panel) and along I = (right panel). Fluxes are averaged over 1° x 1° (dot-dashed, 
orange line), 2° x 2° (continue, red line), 5° x 5° (dashed, brown line) angular bins. 
The corresponding 7-ray flux can be obtained by multiplying this diagram by 3.1. 
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In figure M we compare our results with those obtained in [21] and [22J which have 
been derived assuming a uniform CR density. We also show the flux as obtained by 
using the same gas model B but a uniform CR density. By comparing those profiles the 
reader can see as the flux that we expect from the central region of the GP is significantly 
larger than in [21] . That discrepancy is due, for a large fraction, to the SNR overdensity 
in the molecular ring regiorjjj]. 
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Figure 9. The + flux for E > 1 TeV, averaged over 1° x 1° angular bins, is 
shown. The continuos (red), dot-dashed (blue), dashed (green), and dotted (orange) 
curves correspond respectively to: our model 3B, a model with the same gas density 
distribution (model B) but a uniform CR density as observed at the Earth position; 
the model considered in |21j ; the model considered in [22] . The corresponding 7-ray 
flux can be obtained by multiplying this diagram by 3.1. 

Finally in figure [TOl we represent a full sky map of § u ^ + p^(E u > 1 TeV; bA) obtained 
using model 3B. This figure has been done with the HEALPix package |62jj| 

We conclude this section by discussing how our predictions would change if different 
models for the gas and CR distributions are adopted. We find that by using the gas 
model A rather than B, the 7-ray and neutrino fluxes along the GP grow at most by a 
factor 2. We find comparable displacements by adopting any other of the CR models 
considered in section [3J Therefore, we conclude that, under the hypothesis on which 
this analysis is based, our results are very little model dependent. 



5. Comparison with available observations 

5.1. Gamma-rays 

First of all, we would like to compare our results with EGRET observations. That is 
required to put our results onto more solid grounds and to verify that our choice of the 

|| The reader should also note that in [21] the Nh — » > 1 TeV; b, I) conversion factor is 

6.02 x IP" 33 while we have 8.87 x 10" 33 (TeV s sr) -1 . 
If See jhttp: / /healpix.jpl.nasa.gov| 
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Figure 10. The all-sky map of & v +p (E v > 1 TeV) is shown in galactic coordinates 
(the GC is in the center). This map has been obtained with HEALPix [55]. The flux 
distribution is projected in such a way that all pixels in this map correspond to the 
same solid angle. The corresponding 7-ray flux can be obtained by multiplying this 
diagram by 3.1. 

preferred model is consistent with those observations. Here we consider only EGRET 
measurements of the diffuse emission along the GP between 4 and 10 GeV [5], [63] . Such 
a comparison requires a considerable extrapolation of our previous findings. We are 
allowed to do that since, already for < 10 TeV, nuclei propagation takes place deep into 
the spatial diffusion regime so that the energy dependence of the diffusion coefficients 
do not change going to lower energies (see e.g. [261 I2Z])- As a consequence, in the 
stationary regime the CR distribution above 10 GeV is given by a uniform oc E~ 2 ' 7 
rescaling of that shown in figure [3] We neglect small corrections due to the variations 
of the pp scattering cross section with the energy. 

In table [2] and in figure (TT] we compare respectively theoretical mean fluxes and 
flux profiles along the GP with EGRET measurements [63]. Considering the large 
uncertainties involved in the modelling of the CR and gas distribution, the agreement 
is rather good. It should be noticed that in the energy range 4 — 10 GeV the IC 
contribution to the photon flux from the GP is expected to be subdominat (see e.g. 
[511 |24j). In figure [121 we also show our prediction for the 7-ray differential spectrum 
between 1 GeV and 100 TeV in the region 73.5° < I < 76.5°, < 1.5° and we compared 
it with EGRET [63] and MILAGRO measurments jllj in the same region 0. All the 

+ The best fit of the spectral slope observed by EGRET above 1 GeV is —2.6 [4]. That is compatible 
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Figure 11. The profile of 7-ray flux along the Galactic Plane, between 4 and 10 
TeV, as obtained with our reference model 3B ( blue, dashed line) is compared with 
that measured by EGRET (red, continuos line). Fluxes are averaged over 0.5° x 0.5° 
angular bins. 

theoretical fluxes shown in this section have been derived using the CR model 3 and gas 
model B. As we mentioned in the previous section, while we prefer gas model B because 
it provides the best description of CO surveys, once made that choice, models 3 is that 
which gives the best fit of EGRET data. Anyhow, passing from one model to another, 
our predictions change by a factor 2 at most. 

Reassured by these findings we are now ready to compare our results with 
measurements performed at higher energies. Several air shower array experiments looked 
for the 7-ray diffuse emission of the Galaxy above the TeV. The most interesting results 
are those of TIBET [10] and MILAGRO [TT]. Both experiments probed different 
regions of the GP. In all those regions only upper limits were found but in Cygnus 
where MILAGRO found a significant excess on the background [9j [TT]. In table [2] we 
compare our predictions, as obtained with our preferred 3B model (see section @J, with 
those measurements. With the exception of Cygnus (see also the discussion about the 
GC ridge at the the of this subsection) in all other regions we predict fluxes which are 
significantly below the experimental limits. Therefore, there is still room for a large IC 
contribution. 

Concerning the excess in the Cygnus region observed by MILAGRO, we confirm 
the conclusion [TT], [63] that it cannot be explained by the interaction of the diffuse 
component of galactic CR with the gas in that region (see figure IT2"j) . Indeed, we find 
that a CR local over-density of about 20 is required to explain that signal in terms of 
hadronic emission. 

with that adopted in our paper within errors. 
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Figure 12. The blue (continuos line) is the differential 7-ray spectrum multiplied by 
E 2 as obtained with our reference model 3B in the region 73.5° < I < 76.5°, \b\ < 1.5° 
. Green (red) bars represent EGRET (MILAGRO) measurements in the same window. 

Table 2. In this table our predictions for the mean 7-ray flux in some selected regions 
of the sky are compared with some available measurements. Since measurements' 
errors are much smaller than theoretical uncertainties they are not reported here. 



sky window 




$ 7 (> Ej) ( 


cm 2 s sr) 1 






our model 


measurements 


\l\ < 10°, \b\ < 2° 


4 GeV 


~ 4.7 x 1CT 6 


- 6.5 x 10" 6 [53] 


20° < I < 55°, \b\ < 2° 


3 TeV 

4 GeV 


~ 5.7 x 1CT 11 
~ 4.4 x 1CT 6 


< 3 x 10" 10 [ID] 
- 5.3 x 10" 6 [B3J 


73.5° < I < 76.5°, \b\ < 1.5° 


12 TeV 
4 GeV 


~ 2.9 x 10" 12 1 
~2.4xl0- 6 1 


± 6.0 x 10" 11 [TT] 
t 3.96 x 10" 6 [S3] 


140° < I < 200°, |6| < 5° 


3.5 TeV 
4 GeV 


~ 5.9 x 10" 12 
~ 5.9 x 1CT 7 


< 4 x 1CT 11 [9] 
- 1.2 x 1CT 6 [63] 



A similar effect has to be invoked to explain HESS measurement of the diffuse 
7-ray TeV emission from the GC ridge region [65]. In fact, even accounting for the 
detailed 3-dimensional gas model recently developed in [16] for the GC region, we find 
that in the region |fe| < 0.8°, |/| < 0.3° the expected 7-ray flux for E > 1 TeV is 
~ 2.4 x 10~ 9 cm -2 s _1 sr -1 , which is almost 10 times smaller than HESS measurement. 
It is worth to notice that the slope of the spectrum measured by HESS is close to —2.3 
which is quite different from that of galactic CR. That probably means that primary 
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particles in the GCR have a local origin. It is evident that the analysis performed in 
this work, as well as in related papers [23j EU [21], cannot account for those kind of local 
emissions. 

5.2. Neutrinos 

The only available upper limit on the neutrino flux from the Galaxy has been obtained by 
the AMANDA-II experiment pf] . Being located at the South Pole, AMANDA cannot 
probe the emission from the GC. In the region 33° < I < 213°, \b\ < 2°, and assuming 
a spectral index a = 2.7, their present constraint is 

< 6.6 x 1(T 4 ( n | V ) " (GeV cm 2 s sr)" 1 , (17) 

which implies 1 TeV) < 3.1 x 10~ 9 (cm 2 s sr) -1 . According to our model the 

expected flux in the same region is $i/ M +p M (> 1 TeV) ~ 4.2 x 10~ n (cm 2 s sr) -1 . 

6. The expected neutrino signal in the North Hemisphere 

In the energy range 1 — 100 TeV water (or ice) Cherenkov neutrino telescopes are 
best suited to look for up-going muons produced by charged-current interactions of 
muon neutrinos in the Earth. While the Earth offers almost a complete shielding 
from up-going atmospheric muons below the horizon, an unavoidable background is 
given by atmospheric neutrinos. Several estimates have been made of the atmospheric 
neutrino flux, based on different assumptions modelling hadronic interactions (see e.g. 
[661 [671 [68] ) . Above the TeV all calculations almost agree predicting an averaged flux 

F^{E V ) ~ 4.6 x 10- 8 (r^y) '''TeV -1 cm" 2 s" 1 sr" 1 (18) 

though a ~ 40 % uncertainty remains due to the experimental error on the primary CR 
spectrum and the theoretical error modelling strange particle production. 

Since, as it is evident from our previous results, the expected neutrino flux from 
the Galaxy is significantly smaller than such background, a suitable procedure has to be 
adopted to disentangle the signal. One possible approach is to search only for neutrinos 
with E > 100 TeV [25J. That however may be hard to do due to the very low flux at that 
high energies and to the Earth opacity (only Earth skimming neutrinos or shower-like 
events produced by down-going neutrinos can be detected in that case). In our opinion 
a more promising strategy is to search for up-going muon neutrinos above 1 — 10 TeV. 
Their arrival direction can be reliably reconstructed with an angular resolution as good 
as 0.5° in water and that information may be used to identify the galactic emission as 
a localised excess of events. 

While it is quite evident that the most natural search window (on-source region) 
is a narrow band along the galactic equator (see figure [TO]) . its optimal sizes have to 
be chosen by taking into account the angular extension of the source. In the previous 
section we showed that most of the neutrino flux from the Galaxy should be concentrated 
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in the region \l\ < 50°, \b\ < 1°. Since the Gaussian width of the signal a S i g ~ 2° is 
comparable to the expected angular resolution of neutrino telescopes, some attention 
has to be paid when choosing the latitude width of the search window. By assuming 
that both the line spread function of the experiment and the signal profile along b are 
Gaussian having widths o~i s f and a S i g , respectively, the optimal search window width is 

approximatively given by (see e.g. HU) Ab ~ (1 - 2) (af sf + a*) . Since for a km 3 



'Isf ^ "sigj 

water NT the expected line spread function is o~i s f — 0.5°, we think that Ab = 3° is a 
reasonable value to adopt. 

Although to perform a detailed calculation of the expected signal in a given NT is 
beyond the purposes of this work, here we perform a simplified estimate which is only 
intended to give the reader a feeling of the chances that forthcoming experiments have 
to achieve a positive detection. The expected muon detection rate of neutrinos coming 
from that window is 

Nv{>E v )=f [ [ dEA° s (E)v(b,l) dN » { f; lM (19) 

JAl JAb JE V dE v 

where A eS (E v ) is the effective area of the neutrino telescope and v(l,b) is the visibility 
function (i.e. the fraction of time that a point in a sky with galactic coordinates 
(l,b) spends above the visibility horizon). For example, v(0, 0) = 0.67 at ANTARES 
geographical position. 

Since the effective area is relatively weakly dependent on the arrival direction of 
the neutrinos, we adopt here a mean value obtained by averaging over all possible nadir 
angles. As a reference we use A cS (E l/ ) as provided by the ANTARES collaboration 
(see e.g. [69j) and assume that for a km 3 experiment it will be 40 times larger. This is 
a rapidly growing function of the energy for E < 100 TeV. It is interesting to observe 

that the product A eS (E) — ^ ' — - — , which is related to the detection efficiency, is 

dE 

peaked at about 1 TeV for a a = 2.7 power law spectrum. 

In figure [13] we show the integrated muon neutrino detection rate as a function of the 
minimum energy to be expected in a km 3 NT placed at the same geographical position of 
ANTARES. The upper curve delimiting the signal band corresponds to model 3A, while 
the lower to model 2B. We compare that rate with the atmospheric neutrino detection 
rate expected in ANTARES [70] multiplied by a factor 40 to account for the larger 
effective area. Our calculation accounts for the angular dependence of the expected 
atmospheric neutrino detection rate. 

According to figure [13], the background is about 50 times larger than the signal 
above 1 TeV in the \l\ < 30° region, and about 5 times above 10 TeV. Therefore, an 
excess may be detectable for E > 10 TeV only after a considerable number of years. 

* Here we use an effective area as provided by the ANTARES collaboration for point-like sources. 
According to a recent analysis [70], which adopts a suitable rejection strategy for the residual 
atmospheric muon background, the effective area for an all-sky diffuse flux differs very little from 
that for point-like sources. A dedicated analysis for the kind of emission considered in this work has 
not been performed yet. 
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Figure 13. The expected + detection rate (red band, between continue* 
line) from the regions \l\ < 50°, |6| < 1.5° (left), \l\ < 30°, |6| < 1.5° (centre), 
|Z| < 10°, \b\ < 1.5° (right), is compared with the background due to atmospheric 
neutrinos from the same region (ligth blue band, between dashed lines). An effective 
area 40 times larger than ANTARES has been adopted. E is the neutrino energy. 

7. Conclusions 

We have performed a comprehensive calculation of the neutrino and 7-ray flux 
distributions which are expected to be originated by the interaction of CR nuclei with 
the interstellar gas. 

Our computation improves and updates previous analyses under several aspects. 
Modelling the spatial distribution of CR in the Galaxy, we accounted for a spatial 
dependence of the diffusion coefficients. That approach allowed us to test how different 
models of the turbulent component of the GMF affect the CR distribution and the 
secondary 7-ray flux and to choose our preferred model as that which best reproduces 
EGRET 7-ray observations above few GeVs. Interestingly, we found that such a model 
gives also the best agreement with the hard 7-ray spectra observed for several SNRs. 
Having assumed that CR sources are distributed like SNRs, we adopted a detailed model 
of the spatial distribution of those objects. 

Concerning the gas distribution, we considered two models, both for the HI and H 2 , 
which we applied for the first time in this framework. Although these models have 
significant differences in some regions, we showed that they give rise to neutrino and 7- 
ray fluxes which, when averaged over windows larger than few degrees squared, coincide 
within a factor of a 2. We also found that the column density profile across the GP is 
almost the same for those two models, while it differs significantly from that adopted in 
previous works [211 EE] • As a consequence, we predict the neutrino and 7-ray emissions 
to be more narrowly peaked along the GP. This improves significantly the perspectives to 
disentangle the signal from the almost isotropic background if the experimental angular 
resolution will be better than 1° — 2°. 

We compared our predictions with the available observations. 
Concerning the 7-rays, we found that above the TeV our predicted fluxes are below the 
experimental limits by less than an order of magnitude. Since our analysis does not 
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account for a possible IC contribution, which may be not negligible above the TeV, it 
is possible that those experiments will get a positive signal in the next few years. The 
comparison of the experimental results with our predictions, as well as with those of 
previous works (see e.g. [23J, EH [21]), will help to understand which fraction of that 
emission is of leptonic origin. 

Concerning the signal observed by MILAGRO in the Cygnus region [9j [TTJ [64], we 
confirm that it cannot be explained without invoking a significant CR overdensity, or a 
large IC contribution, which may be due to a local concentration of sources. A similar 
effect takes place in the GC region. We found that in order to explain the excess observed 
by HESS from the GC ridge [65] the CR flux must exceed our prediction by a factor of 
about 10. This is not unexpected since the distribution of star forming regions, hence of 
CR sources, is known to be quite clumpy. It is understood that the kind of simulation 
performed in this work, as well as in previous ones [23l l5T| 124"] , can only model the mean 
CR distribution smoothed on scale of several hundred parsecs. 

We also extrapolated our results down to few GeVs, where the IC scattering 
emission is expected to be subdominant. We found that in the most dense regions 
the flux distribution agrees reasonably well with EGRET observations. This result 
encourages us to improve the accuracy of our analysis by using, for example, more 
detailed models of the gas distribution, so to be able to simulate what GLAST may 
observe above the GeV at least in some limited regions of the sky. 

Going back to neutrinos, we compared our predictions for the muon neutrino flux 
from the GP with the experimental limit recently established by AMANDA-II [14] . 
We found that our predicted flux is almost two orders of magnitude below that limit. 
Unfortunately, also ICECUBE [15] will hardly be able to get a positive signal. 
Since a neutrino telescope placed in the North hemisphere may have better changes, we 
investigated this possibility in some details. Assuming that such an instrument will be 
placed at the same position of ANTARES [16] and have a 40 times larger effective area 
we estimated the expected signal and the background along the Galactic Equator. We 
found that the detection of the smooth component of the diffuse emission may require 
more than 10 years of data taking. It should be noted, however, that our analysis 
provides only a lower limit to the expected emission from the Galactic plane. As we 
mentioned in the above, several observations suggest that the CR and gas distributions 
may be more clumpy than what considered in this work. Furthermore, fluctuations 
of those quantities are likely to be spatially correlated. This may lead to a significant 
enhancement of the neutrino flux from some regions as may be the case for the GC ridge 
[56| 154] and the Cygnus region [TIJ E2] • In a forthcoming paper we will try to develop 
our analysis in order to model at least some of these effects. 
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Appendix A. Models of the gas distributions 

The HI and H2 density distributions are generally inferred from radio observations of 
the 21 cm line and 12 CO rotational line emissions respectively. The galacto-centric 
radial density profiles are estimated by converting the line-of-sight velocity (which is 
determined from the Doppler shift) into heliocentric distance by means of the Galactic 
rotation curve (see e.g. [32]). This operation is impossible for the whole Galaxy and 
it always involves large uncertainties so that some phenomenological models have to be 
invoked. 

We are going here to describe in detail the models sketched in section 12.31 

Model A It has been developed by Nakanishi and Sofue (NS) in [44J for the HI and 
by the same authors in [45J for the H 2 . For both gas components NS developed 
3-D models (r, z, <fi) which, however, do not cover the entire Galaxy. By averaging 
over the azimuth they get 2D-distributions which are approximatively symmetric 
with respect to the galactic plane. NS assumed that the vertical profiles of the 
HI and H2 densities have the form n(r, z) = n(r, 0) sech 2 {log(l + V2) z/z 1/2 (r)] 
which is almost coincident with a Gaussian peaked at the galactic plane half- 
width h(r) = Zi/2(r)/ ^ln(2). The quantity 2 • Zi/2(r) is usually called the Full 
Width to Half Maximum (FWHM) scale height. Then, they derived binned radial 
distributions for zy^{r) and the gas density n(r, 0) along the Galactic Plane. 

Model B It is the product of the combination of results of different analyses which 
have been separately performed for the disk and the galactic bulge. For the H 2 
and HI distributions in the bulge we use a detailed 3D model recently developed by 
Ferriere et al. [16] on the basis of several observations. In that model both the H 2 
and the HI are concentrated in two main structures. The central molecular zone 
(CMZ) appears as a quite dense 500 x 30 pc wide (sizes at half maximum density) 
ellipse on the plane of the sky (for the HI the vertical extension is 90 pc) containing 
almost 2 x 10 7 M Q in H 2 and 1 x 10 6 M Q in HI. It gives rise to a pronounced peak 
in the gas column density in the direction of the GC. The holed Galactic bulge disk 
(GBD) is 3 kpc long and 1 kpc wide toroidal structure. Its mass is comparable to 
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that in the CMZ (~ 3 x 10 7 M in H 2 and 3 x 10 6 M in HI) but it is spread over a 
much larger volume so that the gas density in that region is much smaller than in 
the CMZ. Furthermore, the GBD is significantly inclined with respect to the plane 
of the sky so that its contribution to the gas column density along the line of sight 
is quite small. 

For the molecular hydrogen in the disk we use the Bronfman's et al. model |47j . 
Although such a model is almost two decades old, it was shown [32] that it still 
provides a good description of recent CO surveys [49J. Once averaged over the 
azimuth, Bronfman's et al. distribution is well approximated by a disk with a 
Gaussian vertical profile symmetric with respect to the Galactic plane. Since in 
[4"T] r = 10 kpc was adopted, we correct the gas densities and the scale heights 
given in that paper to make them compatible with the value r = 8.5 kpc we use 
in this work. 

Furthermore, we accounted for the different values of the H 2 -CO conversion factor 
which is X = 0.5 x 10 20 cm~ 2 K" 1 km _1 s in pE] and X = 2.8 x 10 20 cm" 2 K" 1 km^s 
in [17] . Even by taking into account a possible increasing of X with r (see e.g. 
eq.(7) in [46J), these values are too different to be compatible. We assume that 
for 2 < r < 10 X = 1.2 x 10 20 cm -2 K _1 km _1 s as this is the mean value which 
one gets in that region by taking X = 0.5 at r = and assuming that it grows 
like exp(r/7.1) as argued in [73]. Hence, we correct Bronfman's et al. density by 
multiplying it by the factor (1.2/2.8). 



Appendix B. Numerical solution of the diffusion equation 

We will derive here equation ( fill . 

The diffusion equation which describes the propagation of CRs in a turbulent 
medium is ( ITU1) . 

ViJ^E^z) = -V t (D^r,z)VjN(E,r,z)) =Q(E,r,z). (B.l) 

Cylindrical symmetry will be assumed. 

The diffusion tensor can be conveniently decomposed into 

Da = (Dl - D^bibj + D||<J y + D A e ijk b k (B.2) 

where bi are the components of the regular magnetic field versor. The symmetric 
components D\\ and D± are the diffusion coefficients along and perpendicularly to B reg , 
while Da is the antisymmetric (Hall) diffusion coefficient. All those coefficients are 
functions of the energy though with different behaviours. For a fixed value of L max 
and for <C L max , D\\ and D± are proportional to E 2 ^ 1 while Da oc E (see e.g. 
[50. 53J), where 7 is the power-law index of the GMF turbulent fluctuations as defined 
in section 12.21 The Hall diffusion becomes dominant only for very high values of the 
rigidity, i.e. for r L (E) > 0.1L max . 

As motivated in section [2721 and similarly to what done in other works J50J [23j [27J , we 
can approximate the regular component of the magnetic field to be azimutally oriented, 



Diffuse Neutrino and Gamma-ray Emissions of the Galaxy above the Te V 



26 



i.e. = ±1 and b r = b z = 0. Under such an assumption Dm becomes not physically 
relevant and the diffusion equation takes the simpler form [50] 



--d r [rD ± d r ] - d z [D ± d z ] + u r d r + u z d z ^ N(E, r, z) = Q(E, r, z) , (B.3) 



where we defined 

u r = -d z D A (B.4) 

u z = -d r (D A r) . 
r 
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